require(forecast)
require(runjags)
require(coda)
require(stringi)

source("A10_functions.r")

data <- read.csv("C1_main_data.csv")
data$date <- as.Date(data$date)

# augmented Dickey-Fuller test for the LDP support rate
ndiffs(data$LDP, test = "adf")
# Kwiatkowski-Phillips-Schmidt-Shin test for the LDP support rate
ndiffs(data$LDP, test = "kpss")

#### trend in LDP public support rate ####
# Figure 1
png("Figure_1.png", 
    width = 6, height = 3, units = "in", pointsize = 7, res = 1500)
par(mar = c(4, 4, 1, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(1, nrow(data)), ylim = c(0, 40), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = seq(0, 40, 10), lty = 3, col = "gray80")
lines(1:nrow(data), data$LDP)
mtext("LDP support ratings", 2, 2.5)
axis(1, at = c(1, seq(7, 199, 12), nrow(data)), 
     labels = c("Jul 00", "Jan 01", "Jan 02", "Jan 03", "Jan 04", "Jan 05", 
                "Jan 06", "Jan 07", "Jan 08", "Jan 09", "Jan 10", "Jan 11", 
                "Jan 12", "Jan 13", "Jan 14", "Jan 15", "Jan 16", "Jan 17", 
                "Oct 17"), 
     lwd = 0.5, las = 2)
axis(2, lwd = 0.5)
dev.off()

#### estimation ####
# preprocessing
n <- nrow(data) - 1
m <- 4
y <- cbind(data$K_L_3s, data$K_D_3s, data$K_C_3s, data$K_K_3s)[-1, ]
k <- rowSums(y)
y[k == 0, ] <- NA
k[k == 0] <- 1

x <- cbind(diff(data$cabinet / 100), 
           diff(data$LDP / 100), 
           data$DPJ.gov[-1], 
           diff(data$cabinet / 100) * data$DPJ.gov[-1], 
           data$HOR[-1] / 47, 
           (data$HOR[-1] / 47) ^ 2, 
           data$HOC[-1] / 35, 
           (data$HOC[-1] / 35) ^ 2, 
           diff(data$pmid))
l <- ncol(x)

cabinet.mean.minus.SD <- mean(x[x[, 9] == 0, 1]) - sd(x[x[, 9] == 0, 1])
cabinet.mean.plus.SD <- mean(x[x[, 9] == 0, 1]) + sd(x[x[, 9] == 0, 1])
LDP.mean.minus.SD <- mean(x[, 2]) - sd(x[, 2])
LDP.mean.plus.SD <- mean(x[, 2]) + sd(x[, 2])

# initial values
initial.values <- list()
for (i in 1:3) {
  initial.values[[i]] <- inits.fun(10 * i, 100 * i, i)
}

# estimation
result <- run.jags(model = "B1_dynamic_multinomial_model.R", 
                   monitor = c("beta", "Omega", "theta", "theta.star"), 
                   data = list(y = y, k = k, x = x, n = n, m = m, l = l, 
                               mu0 = rep(0, m - 1), 
                               Omega.star0 = diag(m - 1) * 0.01), 
                   inits = initial.values, modules = "glm", 
                   n.chains = 3, burnin = 5000, sample = 5000, 
                   adapt = 20000, summarise = FALSE, 
                   thin = 50, method = "parallel")
# save(result, file = "D1_Komeito_main_result.Rdata")
# load("D1_Komeito_main_result.Rdata")

# convergence check
sum(gelman.diag(result$mcmc, multivariate = FALSE)$psrf[, 1] > 1.1)

# combine MCMC draws into a single matrix
sims.matrix <- rbind(result$mcmc[[1]], result$mcmc[[2]], result$mcmc[[3]])
colnames(sims.matrix) <- colnames(result$mcmc[[1]])

#### results ####
var.names <- c("Change in cabinet approval", "Change in LDP support", 
               "DPJ government", "Cabinet x DPJ gov.", 
               "Time since HoR election", 
               "Time since HoR election squared", "Time since HoC election", 
               "Time since HoC election squared", "Inauguration of the PM")

# Table A.3
round(posterior.summary(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], 
                        var.names), 3)
round(posterior.summary(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], 
                        var.names), 3)
round(posterior.summary(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], 
                        var.names), 3)

#### simulation for the LDP support rate ####
LDP.seq <- seq(round(min(x[, 2]), 3), round(max(x[, 2]), 3), 0.001)

# Figure 2
sim.result.LDP <- array(NA, c(length(LDP.seq), 4, 3))
for (i in 1:length(LDP.seq)) {
  sim.x <- x
  sim.x[, 2] <- LDP.seq[i]
  exp.theta.star.1 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),1\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], sim.x))
  exp.theta.star.2 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),2\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], sim.x))
  exp.theta.star.3 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),3\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], sim.x))
  sum.exp.theta.star <- exp.theta.star.1 + exp.theta.star.2 + exp.theta.star.3 + 1
  prob.1 <- exp.theta.star.1 / sum.exp.theta.star
  prob.2 <- exp.theta.star.2 / sum.exp.theta.star
  prob.3 <- exp.theta.star.3 / sum.exp.theta.star
  prob.4 <- 1 / sum.exp.theta.star
  sim.result.LDP[i, 1, ] <- posterior.summary(rowMeans(prob.1))
  sim.result.LDP[i, 2, ] <- posterior.summary(rowMeans(prob.2))
  sim.result.LDP[i, 3, ] <- posterior.summary(rowMeans(prob.3))
  sim.result.LDP[i, 4, ] <- posterior.summary(rowMeans(prob.4))
}

png("Figure_2.png", 
    width = 6, height = 1.6, units = "in", pointsize = 7, res = 1500)
layout(matrix(1:4, 1, 4))
par(mar = c(4, 4, 2, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.053, 0.089), ylim = c(0.05, 0.25), 
     main = "LDP", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = LDP.mean.minus.SD, lty = 3, col = "gray80")
abline(v = LDP.mean.plus.SD, lty = 3, col = "gray80")
line.fun(LDP.seq, sim.result.LDP[, 1, ])
rug(x[, 2])
box()
mtext("Change in LDP support", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.04, 0.08, 0.02), 
     labels = c("-0.04", "", "0.00", "", "0.04", "", "0.08"), lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.053, 0.089), ylim = c(0, 0.2), 
     main = "DPJ", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = LDP.mean.minus.SD, lty = 3, col = "gray80")
abline(v = LDP.mean.plus.SD, lty = 3, col = "gray80")
line.fun(LDP.seq, sim.result.LDP[, 2, ])
rug(x[, 2])
box()
mtext("Change in LDP support", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.04, 0.08, 0.02), 
     labels = c("-0.04", "", "0.00", "", "0.04", "", "0.08"), lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.053, 0.089), ylim = c(0, 0.2), 
     main = "JCP", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = LDP.mean.minus.SD, lty = 3, col = "gray80")
abline(v = LDP.mean.plus.SD, lty = 3, col = "gray80")
line.fun(LDP.seq, sim.result.LDP[, 3, ])
rug(x[, 2])
box()
mtext("Change in LDP support", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.04, 0.08, 0.02), 
     labels = c("-0.04", "", "0.00", "", "0.04", "", "0.08"), lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.053, 0.089), ylim = c(0.4, 0.8), 
     main = "Komeito", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = LDP.mean.minus.SD, lty = 3, col = "gray80")
abline(v = LDP.mean.plus.SD, lty = 3, col = "gray80")
line.fun(LDP.seq, sim.result.LDP[, 4, ])
rug(x[, 2])
box()
mtext("Change in LDP support", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.04, 0.08, 0.02), 
     labels = c("-0.04", "", "0.00", "", "0.04", "", "0.08"), lwd = 0.5)
axis(2, lwd = 0.5)
dev.off()

# difference in predicted probabilities between cases with changes in the LDP public support rate one SD above vs. below its mean
sim.diff.LDP <- function(low.value, high.value, prob = 0.95) {
  sim.x <- x
  sim.x[, 2] <- low.value
  exp.theta.star.1 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),1\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], sim.x))
  exp.theta.star.2 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),2\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], sim.x))
  exp.theta.star.3 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),3\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], sim.x))
  sum.exp.theta.star <- exp.theta.star.1 + exp.theta.star.2 + exp.theta.star.3 + 1
  prob.1.low <- exp.theta.star.1 / sum.exp.theta.star
  prob.2.low <- exp.theta.star.2 / sum.exp.theta.star
  prob.3.low <- exp.theta.star.3 / sum.exp.theta.star
  prob.4.low <- 1 / sum.exp.theta.star
  sim.x[, 2] <- high.value
  exp.theta.star.1 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),1\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], sim.x))
  exp.theta.star.2 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),2\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], sim.x))
  exp.theta.star.3 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),3\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], sim.x))
  sum.exp.theta.star <- exp.theta.star.1 + exp.theta.star.2 + exp.theta.star.3 + 1
  prob.1.high <- exp.theta.star.1 / sum.exp.theta.star
  prob.2.high <- exp.theta.star.2 / sum.exp.theta.star
  prob.3.high <- exp.theta.star.3 / sum.exp.theta.star
  prob.4.high <- 1 / sum.exp.theta.star
  prob.1.diff <- prob.1.high - prob.1.low
  prob.2.diff <- prob.2.high - prob.2.low
  prob.3.diff <- prob.3.high - prob.3.low
  prob.4.diff <- prob.4.high - prob.4.low
  sim.diff <- rbind(posterior.summary(rowMeans(prob.1.diff), prob = prob), 
                    posterior.summary(rowMeans(prob.2.diff), prob = prob), 
                    posterior.summary(rowMeans(prob.3.diff), prob = prob), 
                    posterior.summary(rowMeans(prob.4.diff), prob = prob))
  rownames(sim.diff) <- c("LDP", "DPJ", "JCP", "Komeito")
  colnames(sim.diff) <- c("Est", "Upper", "Lower")
  cor.matrix <- cor(cbind(rowMeans(prob.2.diff), 
                          rowMeans(prob.3.diff), 
                          rowMeans(prob.4.diff)))
  rownames(cor.matrix) <- colnames(cor.matrix) <- c("DPJ", "JCP", "Komeito")
  list(sim.diff = sim.diff, cor.matrix = cor.matrix)
}
sim.diff.LDP.out <- sim.diff.LDP(LDP.mean.minus.SD, LDP.mean.plus.SD)
round(sim.diff.LDP.out$sim.diff, 3)

# correlation between the predicted probabilities (Online Appendix H)
round(sim.diff.LDP.out$cor.matrix, 3)

# equivalent simulation combining DPJ and JCP into a single category (Online Appendix H)
sim.diff.LDP.combined <- function(low.value, high.value, prob = 0.95) {
  sim.x <- x
  sim.x[, 2] <- low.value
  exp.theta.star.1 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),1\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], sim.x))
  exp.theta.star.2 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),2\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], sim.x))
  exp.theta.star.3 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),3\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], sim.x))
  sum.exp.theta.star <- exp.theta.star.1 + exp.theta.star.2 + exp.theta.star.3 + 1
  prob.1.low <- exp.theta.star.1 / sum.exp.theta.star
  prob.2.low <- exp.theta.star.2 / sum.exp.theta.star
  prob.3.low <- exp.theta.star.3 / sum.exp.theta.star
  prob.4.low <- 1 / sum.exp.theta.star
  sim.x[, 2] <- high.value
  exp.theta.star.1 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),1\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], sim.x))
  exp.theta.star.2 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),2\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], sim.x))
  exp.theta.star.3 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),3\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], sim.x))
  sum.exp.theta.star <- exp.theta.star.1 + exp.theta.star.2 + exp.theta.star.3 + 1
  prob.1.high <- exp.theta.star.1 / sum.exp.theta.star
  prob.2.high <- exp.theta.star.2 / sum.exp.theta.star
  prob.3.high <- exp.theta.star.3 / sum.exp.theta.star
  prob.4.high <- 1 / sum.exp.theta.star
  prob.1.diff <- prob.1.high - prob.1.low
  prob.2.3.diff <- (prob.2.high + prob.3.high) - (prob.2.low + prob.3.low)
  prob.4.diff <- prob.4.high - prob.4.low
  sim.diff <- rbind(posterior.summary(rowMeans(prob.1.diff), prob = prob),
                    posterior.summary(rowMeans(prob.2.3.diff), prob = prob),
                    posterior.summary(rowMeans(prob.4.diff), prob = prob))
  rownames(sim.diff) <- c("LDP", "DPJ/JCP", "Komeito")
  colnames(sim.diff) <- c("Est", "Upper", "Lower")
  sim.diff
}
round(sim.diff.LDP.combined(LDP.mean.minus.SD, LDP.mean.plus.SD), 3)

#### simulation for cabinet approval rate ####
cabinet.seq <- seq(round(min(x[x[, 9] == 0, 1]), 3), 
                   round(max(x[x[, 9] == 0, 1]), 3), 0.005)

# Figure A.5
sim.result.cabinet <- array(NA, c(length(cabinet.seq), 4, 3))
for (i in 1:length(cabinet.seq)) {
  sim.x <- x
  sim.x[, 1] <- cabinet.seq[i]
  sim.x[, 3] <- sim.x[, 4] <- sim.x[, 9] <- 0
  # exclude weeks from 111 to 149 (DPJ government period)
  exp.theta.star.1 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|10[0-9]|110|1[5-9][0-9]|2[0-1][0-9]),1\\]") &
                                                            (! stri_detect_regex(colnames(sims.matrix), paste0("star\\[(", paste(which(x[, 9] == 1), collapse = "|"), "),1\\]")))]), rep(1, nrow(sim.x))) +
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], sim.x))
  exp.theta.star.2 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|10[0-9]|110|1[5-9][0-9]|2[0-1][0-9]),2\\]") &
                                                            (! stri_detect_regex(colnames(sims.matrix), paste0("star\\[(", paste(which(x[, 9] == 1), collapse = "|"), "),2\\]")))]), rep(1, nrow(sim.x))) +
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], sim.x))
  exp.theta.star.3 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|10[0-9]|110|1[5-9][0-9]|2[0-1][0-9]),3\\]") &
                                                            (! stri_detect_regex(colnames(sims.matrix), paste0("star\\[(", paste(which(x[, 9] == 1), collapse = "|"), "),3\\]")))]), rep(1, nrow(sim.x))) +
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], sim.x))
  sum.exp.theta.star <- exp.theta.star.1 + exp.theta.star.2 + exp.theta.star.3 + 1
  prob.1 <- exp.theta.star.1 / sum.exp.theta.star
  prob.2 <- exp.theta.star.2 / sum.exp.theta.star
  prob.3 <- exp.theta.star.3 / sum.exp.theta.star
  prob.4 <- 1 / sum.exp.theta.star
  sim.result.cabinet[i, 1, ] <- posterior.summary(rowMeans(prob.1))
  sim.result.cabinet[i, 2, ] <- posterior.summary(rowMeans(prob.2))
  sim.result.cabinet[i, 3, ] <- posterior.summary(rowMeans(prob.3))
  sim.result.cabinet[i, 4, ] <- posterior.summary(rowMeans(prob.4))
}

png("Figure_A5.png", 
    width = 6, height = 1.6, units = "in", pointsize = 7, res = 1500)
layout(matrix(1:4, 1, 4))
par(mar = c(4, 4, 2, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", xlim = c(-0.221, 0.136), ylim = c(0.05, 0.25), 
     main = "LDP", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = cabinet.mean.minus.SD, lty = 3, col = "gray80")
abline(v = cabinet.mean.plus.SD, lty = 3, col = "gray80")
line.fun(cabinet.seq, sim.result.cabinet[, 1, ])
rug(x[x[, 3] == 0, 1])
mtext("Change in cabinet approval", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.2, 0.1, 0.05), 
     labels = c("-0.2", "", "-0.1", "", "0.0", "", "0.1"), lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", xlim = c(-0.221, 0.136), ylim = c(0, 0.2), 
     main = "DPJ", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = cabinet.mean.minus.SD, lty = 3, col = "gray80")
abline(v = cabinet.mean.plus.SD, lty = 3, col = "gray80")
line.fun(cabinet.seq, sim.result.cabinet[, 2, ])
rug(x[x[, 3] == 0, 1])
mtext("Change in cabinet approval", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.2, 0.1, 0.05), 
     labels = c("-0.2", "", "-0.1", "", "0.0", "", "0.1"), lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", xlim = c(-0.221, 0.136), ylim = c(0, 0.2), 
     main = "JCP", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = cabinet.mean.minus.SD, lty = 3, col = "gray80")
abline(v = cabinet.mean.plus.SD, lty = 3, col = "gray80")
line.fun(cabinet.seq, sim.result.cabinet[, 3, ])
rug(x[x[, 3] == 0, 1])
mtext("Change in cabinet approval", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.2, 0.1, 0.05), 
     labels = c("-0.2", "", "-0.1", "", "0.0", "", "0.1"), lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", xlim = c(-0.221, 0.136), ylim = c(0.6, 0.8), 
     main = "Komeito", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = cabinet.mean.minus.SD, lty = 3, col = "gray80")
abline(v = cabinet.mean.plus.SD, lty = 3, col = "gray80")
line.fun(cabinet.seq, sim.result.cabinet[, 4, ])
rug(x[x[, 3] == 0, 1])
mtext("Change in cabinet approval", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.2, 0.1, 0.05), 
     labels = c("-0.2", "", "-0.1", "", "0.0", "", "0.1"), lwd = 0.5)
axis(2, lwd = 0.5)
dev.off()

# difference in predicted probabilities between cases with changes in cabinet approval one SD above vs. below its mean
sim.diff.cabinet <- function(low.value, high.value, prob = 0.95) {
  sim.x <- x
  sim.x[, 1] <- low.value
  sim.x[, 3] <- sim.x[, 4] <- sim.x[, 9] <- 0
  exp.theta.star.1 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|10[0-9]|110|1[5-9][0-9]|2[0-1][0-9]),1\\]") &
                                                            (! stri_detect_regex(colnames(sims.matrix), paste0("star\\[(", paste(which(x[, 9] == 1), collapse = "|"), "),1\\]")))]), rep(1, nrow(sim.x))) +
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], sim.x))
  exp.theta.star.2 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|10[0-9]|110|1[5-9][0-9]|2[0-1][0-9]),2\\]") &
                                                            (! stri_detect_regex(colnames(sims.matrix), paste0("star\\[(", paste(which(x[, 9] == 1), collapse = "|"), "),2\\]")))]), rep(1, nrow(sim.x))) +
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], sim.x))
  exp.theta.star.3 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|10[0-9]|110|1[5-9][0-9]|2[0-1][0-9]),3\\]") &
                                                            (! stri_detect_regex(colnames(sims.matrix), paste0("star\\[(", paste(which(x[, 9] == 1), collapse = "|"), "),3\\]")))]), rep(1, nrow(sim.x))) +
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], sim.x))
  sum.exp.theta.star <- exp.theta.star.1 + exp.theta.star.2 + exp.theta.star.3 + 1
  prob.1.low <- exp.theta.star.1 / sum.exp.theta.star
  prob.2.low <- exp.theta.star.2 / sum.exp.theta.star
  prob.3.low <- exp.theta.star.3 / sum.exp.theta.star
  prob.4.low <- 1 / sum.exp.theta.star
  sim.x[, 1] <- high.value
  exp.theta.star.1 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|10[0-9]|110|1[5-9][0-9]|2[0-1][0-9]),1\\]") &
                                                            (! stri_detect_regex(colnames(sims.matrix), paste0("star\\[(", paste(which(x[, 9] == 1), collapse = "|"), "),1\\]")))]), rep(1, nrow(sim.x))) +
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], sim.x))
  exp.theta.star.2 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|10[0-9]|110|1[5-9][0-9]|2[0-1][0-9]),2\\]") &
                                                            (! stri_detect_regex(colnames(sims.matrix), paste0("star\\[(", paste(which(x[, 9] == 1), collapse = "|"), "),2\\]")))]), rep(1, nrow(sim.x))) +
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], sim.x))
  exp.theta.star.3 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|10[0-9]|110|1[5-9][0-9]|2[0-1][0-9]),3\\]") &
                                                            (! stri_detect_regex(colnames(sims.matrix), paste0("star\\[(", paste(which(x[, 9] == 1), collapse = "|"), "),3\\]")))]), rep(1, nrow(sim.x))) +
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], sim.x))
  sum.exp.theta.star <- exp.theta.star.1 + exp.theta.star.2 + exp.theta.star.3 + 1
  prob.1.high <- exp.theta.star.1 / sum.exp.theta.star
  prob.2.high <- exp.theta.star.2 / sum.exp.theta.star
  prob.3.high <- exp.theta.star.3 / sum.exp.theta.star
  prob.4.high <- 1 / sum.exp.theta.star
  prob.1.diff <- prob.1.high - prob.1.low
  prob.2.diff <- prob.2.high - prob.2.low
  prob.3.diff <- prob.3.high - prob.3.low
  prob.4.diff <- prob.4.high - prob.4.low
  sim.diff <- rbind(posterior.summary(rowMeans(prob.1.diff), prob = prob), 
                    posterior.summary(rowMeans(prob.2.diff), prob = prob), 
                    posterior.summary(rowMeans(prob.3.diff), prob = prob), 
                    posterior.summary(rowMeans(prob.4.diff), prob = prob))
  rownames(sim.diff) <- c("LDP", "DPJ", "JCP", "Komeito")
  colnames(sim.diff) <- c("Est", "Upper", "Lower")
  sim.diff
}
round(sim.diff.cabinet(cabinet.mean.minus.SD, cabinet.mean.plus.SD), 3)